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We are presenting an Internal Linear Combination (ILC) CMB map, in which the foreground 
is reduced through harmonic variance minimization. We have derived our method by converting a 
general form of pixel-space approach into spherical harmonic space, maintaining full correspondence. 
By working in spherical harmonic space, spatial variability of linear weights is incorporated in a self- 
contained manner and our linear weights are continuous functions of position over the entire sky. 
The full correspondence to pixel-space approach enables straightforward physical interpretation on 
our approach. In variance minimization of a linear combination map, the existence of a cross term 
between residual foregrounds and CMB makes the linear combination of minimum variance differ 
from that of minimum foreground. We have developed an iterative foreground reduction method, 
where perturbative correction is made for the cross term. Our CMB map derived from the WMAP 
data is in better agreement with the WMAP best-fit ACDM model than the WMAP team's Internal 
Linear Combination map. We find that our method's capacity to clean foreground is limited by the 
availability of enough spherical harmonic coefficients of good Signal-to- Noise Ratio (SNR). 

PACS numbers: 98.70.Vc, 98.80.Es 
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I. INTRODUCTION 



A whole-sky map contains significant amount of fore- 
ground emission from astrophysical sources. Hence, the 
ability to clean foreground contamination in CMB data 
is of the utmost importance for CMB observations. In 
the WMAP observation, foreground were cleaned by two 
different methods P, [3, Q. One is using external tem- 
plates of the foregrounds, but using external template 
maps [l[ suffers from dubious extrapolation and template 
noise of higher level than the WMAP Q . The other is us- 
ing the Internal Linear Combination (ILC) method 0, Q , 
where the linear weight for each frequency channel is cho- 
sen to minimize the variance of the linear combination 
of multi- frequency maps, therefore minimizing residual 
foreground. To take into account the spatial variability, 
the WMAP team defined twelve disjoint regions, where 
distinct linear weights are assumed for each region. In 
spite of many merits of the ILC method, it has important 
limits: First, the definition of disjoint regions requires ex- 
ternal information and the use of disjoint regions brings 
about discontinuities. Second, there exists a cross term 
between the residual foreground and CMB, which makes 
the variance minimization proceed as to maximize the 
cancellation between the residual foreground and CMB 
3]. For the solution of the first problem, we have car- 
ried out the variance minimization entirely in spheri- 
cal harmonic space, where the spatial variability of lin- 
ear weights can be incorporated in a self-contained and 
seamless manner. For the solution of the second prob- 
lem, we have developed an iterative foreground reduc- 
tion method, where perturbative correction is made for 



the cross term. Simulations confirmed that our iterative 
method reconstructs the CMB with stability and relia- 
bility. We have also applied our method to the WMAP 
data and obtained a foreground-reduced CMB map. 

The outline of this paper is as follows. In Sec. [ill 
we discuss briefly the foreground reduction method with 
multi-frequency maps, minimum variance principle and 
the choice of general form for linear weights. In Sec. IIII) 
we derive equations in spherical harmonic space, whose 
solutions correspond to linear weights of minimum fore- 
ground. We present an iterative foreground reduction 
method in Sec. IIVI and simulation results in Sec. El 
The result of application to the WMAP three year and 
five year data are presented in Sec. |VI] and IVIII respec- 
tively. We discuss computational issues in Sec. I Villi and 
conclude this investigation in Sec. IIXI In appendix [Al 
the equation for minimum foreground is put in matrix 
notation and the solution is presented in the form of ma- 
trix operations. In appendix [Bj we show that the cross 
term between residual foreground and CMB causes the 
suppression on the lowest multipole powers of an Internal 
Linear Combination (ILC) map. In appendix [Cl we make 
a brief comparison by summarizing the advantages and 
disadvantages of ILC variants and the template-fitting 
method. 



II. FOREGROUND REDUCTION WITH 
MULTI-FREQUENCY MAPS 

With neglect of pixel noise, a thermodynamic temper- 
ature map at a frequency Vi and pixel x is as follows: 



T(x, = T cmb (x) + T fg (x, Vi), 



(1) 
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where T cm b(x) and Tf g (x, vi) are CMB signal and the 
composite foreground signal respectively. A natural 
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choice for the estimator of the CMB map is a linear com- 
bination of multi-frequency maps, which is as follows: 



Wi(x)T(x,z^). 



To keep the CMB unchanged, a contraint is given such 
that the sum of linear weights over frequency channels is 
equal to unity: 



5^u>i(x) = 1. 



(2) 



With Eq. [T] and [2J it is straightforward to show that 
^ w ( (x) T(x, Vi) = T cmb (x) + ^ Wi(-x.) T fg (x, Uj). (3) 

i i 

We can make the foreground signal in Eq. [3] vanish, if 
the linear weights are chosen such that : 



^Wi(x) T fg (x, = 0. 



(4) 



Since we have no information on Tf g (x, Vj,), we need some 
function to maximize or minimize, which will lead us to- 
ward such linear weights. One of such powerful methods 
is variance minimization of the linear combination map 
[H, [f| . It can be shown that the variance of a linear com- 
bination map is 



Wi(x) T(x, v. 




(5) 



« c 2 + 2 ^T cmb (x) ^ u>i(x) r fg ( x , 

+ ^fctfl i (x)T fi! (x,i/ i )j \ 

where the constant term C 2 is the variance of CMB 
and therefore, independent of the choice of linear weight. 
Though the cross term 2 (T cm b(x) ^ i w,(x) 7f g (x, v^)) in 
Eq. [5] vanishes, when averaged over a whole ensemble of 
universes, it is not necessarily zero for our single observ- 
able Universe. Hence, we assume the cross term to be 
small but non-zero, and will make perturbative correc- 
tion for it (see Sec. IIV[) . For now, we neglect the cross 
term. 

The linear weights, which yield a foreground-free map, 
are functions of the frequency spectrum of foreground 
components. Since the frequency spectrum varies over 
sky (see @ for a recent treatment), the linear weights 
should possess spatial variability. To accommodate the 
spatial variability of linear weights, the WMAP team de- 
fined twelve disjoint regions in the WMAP three year 
ILC (WILC3YR) construction, where distinct values of 
linear weights are assumed for each region. The linear 
weights of WILC3YR have the form Wij , where i and j 



denote a frequency channel and a region index. Though 
the WMAP team used regions of smoothed boundaries 
in the final map making, there still exist intrinsic dis- 
continuities from the use of disjoint regions in variance 
minimization, which may even create artificial peculiari- 
ties. 

To reflect the varying powers of foregrounds on differ- 
ent angular scales, linear weights contrived by Tegmark 
et al. has multipole dependency as well [7[, and have 



the form 



We can easily show that optimal linear 



weights should possess m dependency as well as I de- 
pendency. For illustrative purposes, let's consider two 
frequency channel observation and assume the signal to 
consist of CMB and one foreground component only. The 
spherical harmonic coefficient of ith channel is given by 

a \m — a hn h + a im S ' wnere a \m denotes the spherical har- 
monic coefficient of a foreground at ith channel. Keeping 
the CMB signal unchanged, we assign a linear weight w 
and (1 — w) to the frequency channel 1 and 2 respectively. 
Then, the spherical harmonic coefficient of a linear com- 
bination map is given by 



+ (I - w)t 



cmb 



l>fg i n \ 2,fg 

Wa lm + I 1 - W ) a lm ■ 



Obviously the linear weight w yielding a foreground-free 
linear combination map is 



w 



'7 m 



(6) 



As shown in Eq. [SJ optimal linear weights should possess 
m dependency as well as I dependency. 

The linear combination map of minimum foreground 
formed with multi-frequency maps 



(7) 



can be rewritten in the spherical harmonic space, using 
the Clebsch-Gordon relation as: 



a LM 



(8) 



(-1) 



M 



2L + 1 



I I 1 



££v/(2Z + l)(2Z' + l) 

W lm a \'m'i 
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—M j I 



where 



aLM 



Y* M (e, < p)T(e, ( f>)dn, 



Y l r m ,(e, ( f>)T(e, < p,u i )dn. 



The constraint ^2,iW % {9 , <f>) — 1 imposed to preserve the 
CMB signal is expressed in spherical harmonic space as 
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follows: 
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(9) 
(10) 



We can see that linear weights W; m in Eq. [8] possess 
m dependency as well as I dependency. Since Eq. [5] is 
equivalent to [7J physical interpretation on Eq. [5] is quite 
straightforward and we base our approach on Eq. [8j 



III. DETERMINATION OF LINEAR WEIGHTS 

Through variance minimization, we are going to derive 
equations leading toward the linear weights of minimum 
foreground. Since the function w l (9,4>) is real- valued, 



obeys the reality condition w\_ r 



(-1)' 



w 



Therefore, only w\ m (to > 0) needs to be determined. 
It is computationally convenient to accommodate the re- 
ality condition by defining real- valued spherical harmonic 
coefficients wj m as Re[u> ; l m ], Im[w; m ] for m > 0, to < 
respectively. The constraints on wj m derived from Eq. [5] 
and ITD1 are as follows: 



E*oo = v^, 

i 

5>L = o (z>o). 



(11) 

(12) 



The linear weights of minimum foreground minimize 
the variance X^lm l a £M| under the constraints Eq. 
1111 and 1121 The constrained minimization problem is 
solved conveniently via Lagrange's undetermined multi- 
plier method @. With the introduction of Lagrange's 
multiplier A; TO , it can be shown that the variance is min- 
imized under the constraints Eq. [11] and 1121 when 

LM i \ V i / 



9 ™Vm> 



d ™Vm> 



+ E^ 

l>0,m 



0£< 



Im 



I in ' 



0.(13) 



By using Eq. [SI it can be shown that Eq. [T3] has the 
following form: 



E 

Urn 



r i'i 

[ a l> m > 



Im w lm 



A/' m ' — 0, 



(14) 



where ay mHm is 



I'rn'lm 



= 2Rc 



E 

LM 



7*, (/', to', L, M) to, L, M) 



,(15) 



and 7i(/i,mi,Z 3 , to 3 ) is 

ji(h,mi,l 3 , m 3 ) + (-l) mi ji(h, -mi,l 3 , m 3 ) 

7 4 (?l,TOl,Z 3 ,TO 3 ) 

«[7i(Zi,-mi,/ 3 ,TO 3 ) - (-l) mi 7i(/i,TOi,Z 3 ,TO 3 )] 
for mi > 0, mi = and mi < respectively, and 
7i(?i,m,i,Z 3 ,m 3 ) = 



(16) 



12^2 



(2/i + l) (2/ 2 + l) (2/ 3 + l) 



/ 3 



mi TO2 — m 3 



4tt 

/i ^2 h 




Therefore, the values of linear weights of minimum fore- 
ground can be found in terms of Lagrange's multiplier 
Arm' by solving the system of simultaneous linear equa- 
tions given by Eq. [Mj The values of Lagrange's multi- 
plier Xi> m > can be easily determined by making the so- 
lutions of Eq. [14] satisfy the constraints Eq. QT] and [12] 
We can write Eq. [11] [T^] and [T3] in matrix form and ob- 
tain the solution conveniently via matrix operations (For 
details on the solution in matrix notation, refer to Eq. 
A3). 



PERTURB ATIVE CORRECTION FOR THE 
CROSS TERM 



IV. 



There exists a non-zero correlation between fore- 
grounds and true CMB, so called 'Cosmic Covariance' 0], 
which leads to a non-negligible cross term in Eq. [5] The 
existence of this non-negligible cross term makes the lin- 
ear combination of minimum variance differ from that of 
minimum residual foregrounds [H, 0| ■ By noting that the 
cross term disappears in the absence of CMB signal, we 
have developed a perturbative method, where the cross 
term is reduced through iterations. However, as it was 
shown in 0], this approach is only effective down to the 
level of 'Cosmic Covariance' [9|. 

Consider the quantity T(x, Vj) — rf~^(x), where 

^cmb( x ) i s our b es t g uess CMB map from the (j — l)th 
iteration with 



jiO 



crab 



(x) = 0. 



(17) 



The merit of this quantity is that CMB signal is re- 
duced through iterations, leading to reduction of the 
cross term. We obtain linear weights wfix) of the jth 
iteration through variance minimization of 

£X"(x) (T(x,^-Ttb(*)), 



and update our best guess CMB map as follows: 

7tb« = ^w+E^'w ( r (^)-^(| 
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Using Eq. [T]and[2l we may show that our updated CMB 
map is 

^ mb (x) = T cmb (x) + WW T %( x > "0> (18) 

i 

where T cm b(x) is a true CMB map. Using Eq. 
Q] and [2j we may also show that the variance of 

E^'(x) (tm-t^)) is 



= < ^ cmb (x) - T^x) + WW 7> g (x, ^ 
= C 2 + (^W(x)T fg (x,^)^ ) 

+2<(r cmb (x) - X^x)) (j2 WW r £g (x, 



(19) 




where C 2 is a term independent of linear weights W( x )- 
Using Eq. [TH1 the cross term, which is in the last line of 
Eq. 033 may be shown to be 

-2< W^^x) T fg (x, ^ ^ ^(x) T fg (x, ^ ). 

Therefore, the cross term is getting reduced through it- 
erations, provided that 

;W +1 (x)T fg (x,^ < ^W _1 W7f g (x,^ . 

In practice, the cross term converges to some non-zero 
value, which may be attributed to two causes. First, fore- 
ground are reduced with some error, which arises from 
the imperfection of the method applied (e.g. a finite num- 
ber of assumed wj m , imperfection of twelve disjoint re- 
gions of WILC3YR) . The residual foreground related to 
the error are reduced barely through iterations. Second, 
residual foregrounds of the jth iteration possesses some 
level of correlation with that of the j — 1th iteration. 



CMB map. B\ is the beam transfer functions of the 
WMAP jth channel @ and Bj LC is the beam transfer 
function of a 1° FWHM Gaussian beam, which is the 
smoothing kernel used in WILC3YR. Our procedure for 
the generation of simulated data is overly conservative, 
because the presence of instrument noise and residual 
foreground in the WILC3YR makes some foreground and 
instrument noise double-counted. Using Eq. 1201 we have 
generated four hundred simulated data set and carried 
out foreground reduction on them. 





1 r 




0.8 




0.6 




0.4- 




0.2 


o 






0; 


< 


-0.2 




-0.4 




-0.6 




-0.8 




-1 - 



o 



1 

0.8 
0.6h 
0.4 
0.2 


-0.2 
-0.4 
-0.6 
-0.8 
-1 



100 200 300 

realization index 

FIG. 1: 1 = 2 
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FIG. 2: I = 3 
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APPLICATION TO SIMULATED DATA 



We have generated simulated data as follows: 



B. 



"Ira 
i,A 



(20) 



a im ls ^ ne spherical harmonic coefficients of the WMAP 
band maps at ith channel from the LAMBDA site, and 
a Im° an d a Tm are those of WILC3YR and a simulated 



When linear weights are obtained through variance 
minimization on noisy data, linear weights are chosen 
as to minimize noise rather than foreground. Since our 
simulated data are quite noisy on multipoles higher than 
300, we have used only a\ m (I < 300) in variance mini- 
mization (i.e. summation over li was done up to 300 in 

Eq. mi). 

We have made the assumption that linear weights of 
minimum foreground will be spatially coherent on small 
angular scales, and determined only wi m in a finite mul- 
tipole range < I < l cu toS- Though ideally the total 
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merically non-singular, is much smaller than the ideal 
case (i.e. Z cu toff *C 300) (See Section IVIIII for the discus- 
sion on the possible sources of numerical singularity.). 
We have increased i C utoff until the numerical instability 
emerges and found that i C utoff — 7 is optimal for the 
WMAP data. If more a\ of good SNR were available, 
it would be numerically stable with higher Cutoff- 
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FIG. 3: 1 = 4 
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FIG. 6: realization #22 
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We have implemented iterative foreground reduction, 
which is discussed in Sec. IIVI The cross term in the first 
iteration is quite significant, since the best-guess CMB 
map of zeroth iteration is set to zero (see Eq. [P7|) . Hence, 
in the implementation of iterative foreground reduction, 
we have excluded lowest multipoles (0 < I < 10) in vari- 
ance minimization of the first iteration for regularization 
purpose (see appendix [B] for details on why exclusion of 
lowest multipoles reduces the effect of the cross term.). 
Such regularization is not necessary in succeeding itera- 
tions (J > 2). Since it turned out that the improvement 
after j = 2 iterations were negligible, we carried out only 
j = 2 iterations for each simulated data set. 

The power spectra discrepancy between output CMB 
and input CMB are shown for four hundred CMB real- 
izations in Fig. [U [31 g] and [5] on multipoles (2 < I < 
5, I = 20). The vertical axis denotes (Cf ut - Cf )/C ; in , 
where C\ — (21 + 1) _1 J2 m \ a im\ 2 - The horizontal axis 
denotes the enumerating index of four hundred CMB re- 
alizations. In Fig. [SI we show power spectra (2 < I < 30) 
of input CMB and output CMB for the realization #22, 
since the realization #22 among four hundred realiza- 
tions has the octupole and quadrupole power closest to 
those of the WMAP best-fit ACDM. 



FIG. 5: / = 20 



number of w\ m may be as high as the total number of 
available a\ m (i.e. the number of unknowns may be as 
many as the number of constraints), we found that the 
number of w l lm , which keeps the matrices in Eq. IA5I nu- 



VI. APPLICATION TO THE WMAP THREE 
YEAR DATA 

We have applied our foreground reduction method to 
the WMAP three year data from the LAMBDA site 0]. 
Spherical harmonic coefficients of the band maps have 
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been obtained as follows: 




i,A / pi, A 



IB] 



where a 1 ,^ are the spherical harmonic coefficients of the 
WMAP three year band maps and B\ ,A are the beam 
transfer functions of the WMAP ith channel [|[ ■ Just as 
the application to the simulated data in the previous sec- 
tion, the cutoff multipolc Cutoff for linear weights is set 
to seven, and we have used a\ m in the multipole range 
I < 300 in variance minimization (i.e. summation over 
h was done up to 300 in Eq. [TS}. In Fig. our CMB 
map, which we call 'Harmonic Internal Linear Combi- 
nation map' (hereafter, HILC3YR ), is shown with the 
WMAP three year ILC map (WILC3YR) and the differ- 
ence map. The maps in Fig. [7|are images smoothed with 
1° FWHM beam. 
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FIG. 7: the 1° FWHM smoothed maps hiK]: HILC3YR of 
the zeroth iteration (top), HILC3YR of the first iteration (sec- 
ond), WILC3YR (third), WILC3YR - HILC3YR (fourth), 
and the difference between the zeroth and the first iteration 
HILC3YR (bottom) 
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FIG. 8: Power spectra of HILC3YR, WILC3YR and the 
WMAP 3 year best-fit ACDM model 

Power spectra estimate on HILC3YR and WILC3YR 
are made by computing C; = (2Z + 1) -1 J2 m \ a im\ 2 , which 
are shown with the WMAP best-fit ACDM model in 
Fig. H As shown in Fig. |SJ HILC3YR makes better 



agreement with the WMAP best-fit ACDM model than 
WILC3YR (e.g. The first Doppler acoustic peak is vis- 
ible in the HILC3YR power spectrum around I ~ 220). 
The huge excess power of WILC3YR and HILC3YR on 
high multipoles (/ > 300) is attributed to pixel noise. 
In the multipole range (200 < I < 300), where point 
sources are dominant over other sources [lOj, HILC3YR 
makes relatively good agreement with the model, while 
there is significant discrepancy between WILC3YR and 
the model. This may indicates relatively effective reduc- 
tion of point sources in HILC3YR. The absolute value 
of the power spectra difference between WILC3YR and 
HILC3YR is also shown in Fig. [8] Through polynomial- 
fitting in the multipole range (I > 100), we have inves- 
tigated the multipole dependency of the the power spec- 
tra difference between HILC3YR and WILC3YR. Con- 
sidering the multipole dependency, we may, with some 
caution, attribute the power difference to residual point 
sources and pixel noise. We have also investigated the lin- 
ear weights of HILC3YR, and found that HILC3YR gets 
contribution more from the V band map than the W band 
map. Hence, the less noise level of HILC3YR does not 
mean that HILC3YR prefers blindly the least noise chan- 
nel at the sacrifice of foreground reduction. The linear 
weights of HILC3YR, which are continuous over entire 
sky, are shown in Fig. [9l We have computed the vari- 
ance of our linear weights by WJ = (21 + l) -1 J2 m \ w im\ 2 
to quantify the spatial variation of our linear weights on 
different angular scales. As shown in Fig. [TQl Wj tends 
to decrease with increasing multipole, with Wq being the 
highest. It is not difficult to sec from the tail pattern that 
there will be some non-zero w\ m on multipoles higher 
than our assumed cutoff multipole I = 7. These unac- 
counted w\ m (I > 7) may be partially responsible for the 
residual foregrounds, which is visible around the galactic 
plane in Fig. [7] However, it is unlikely that low multi- 
pole anisotropics arc affected significantly by the residual 
foreground around the Galactic plane. 

TABLE I: quadrupole and octupole powers 



Measurement ST^[fiK A ] p- value 8Tg[p,K^] 



WMAP best-fit ACDM model 


1250 




1143 


Hinshaw et al. cut sky 


211.0 


2.6% 


1041 


WMAP team's ILC (WILC3YR) 


248.6 


3.7% 


1051.5 


Tegmark et al. (TCM3YR) 


209.6 


2.5% 


1037.8 


HILC3YR 


331 


6.8% 


961 



In Table HI the quadrupole and octupole power of 
HILC3YR are shown with those of the WMAP best- 
fit ACDM model and other measurements. The p- value 
in Table Q] denotes the chance of having the quadrupole 
power lower than the measurements on the left. As shown 
in Table HI the p-value of HILC3YR quadrupole is al- 
most twice that of WILC3YR. The preferred axis of ar- 
bitrary multipoles can be quantified by finding the axis 
n, which maximizes the angular momentum dispersion 
of the corresponding multipole [ll[ . It was noticed that 
the preferred axis of quadrupole anisotropy is close to 
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FIG. 9: The HILC3YR linear weight for the K, Ka, Q, V, and 
W band map (from top to bottom) 



being in alignment with that of octupole [ll| . The angu- 
lar separation between the preferred axis of quadrupole 
and octupole of HILC3YR is shown with that of other 
foreground-reduced maps in Table HI The p- values de- 
notes the probability of the angular separation lower than 
the measurements, provided that the direction of a pre- 
ferred axis is random. 



TABLE II: quadrupole-octupole alignment 



Maps 


#23 


p-value 


WMAP team's ILC (WILC3YR) 


5° .9 


0.53% 


Tegmark et al. (TCM3YR) 


13°. 2 


2.65% 


HILC3YR 


13°. 1 


2.60% 
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FIG. 10: Variance of the HILC3YR linear weights 

VII. APPLICATION TO THE WMAP FIVE 
YEAR DATA 




FIG. 11: HILC3YR (top), HILC5YRa (middle), HILC3YR - 
HILC5YRa (bottom) 



the preparation of this paper. For comparison with 
HILC3YR, we have obtained HILC5YRa using the same 
HILC parameters as the HILC3YR (i.e. Cutoff = 7 and 
a \m °f the- multipole range (I < 300) for variance mini- 
mization). The HILC5YRa is shown with the HILC3YR 
in Fig. QT] We may see that the HILC5YRa contains less 
level of residual foreground, since most of the difference 
map shown in Fig. [TT] is positive. It is reported that 
some improvement in instrument calibration have been 
made for the WMAP 5 year data [Hj]. We attribute 
less level of residual foregrounds in the HILC5YRa par- 
tially to the instrument calibration improvement, since 
it might improve the accuracy of frequency dependency 
of band map data. We found that Eq. IA5I used with the 
WMAP 5 year data has less degree of numerical singular- 
ity than the 3 year data, which may also be attributed to 
the improved accuracy of frequency dependency of band 
map data. Since the WMAP 5 year data have higher 
SNR than the 3 year data, we have increased the multi- 
pole range of a\ m to I < 400. The improved numerical 
stability and the increase in the multipole range of a\ m 
(I < 400) allowed us to set Z C utoff to 15. Hence we have 
obtained the HILC5YR with the / C utoff = 15, which is 
shown with WILC5YR in Fig. [12] 




We have also applied our method to the WMAP FIG - 12: HILC5YR (top), WILC5YR (middle), WILC5YR - 
five year data [II, which have been released during HILC5YR (bottom) 
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The power spectra estimate on HILC5YR and 
WILC3YR are shown with the WMAP 5 year best-fit 
ACDM model in Fig. [131 By comparing Fig. [13] with 
Fig. [51 we may see that noise level in 5 year ILC maps 
are lower than that of 3 year ILC maps, as expected. It 
was reported that the existence of the cross term leads to 
the suppression of low multipole anisotropy of ILC maps 
[iltljl- It is interesting to note that most of low multipole 
powers of WILC5YR are lower than those of HILC5YR. 
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FIG. 13: Power spectra of HILC5YR, WILC5YR and the 
WMAP best-fit ACDM model 



TABLE III: quadrupole and octupole powers obtained with 
the WMAP 5 year data 



Measurement 




p- value 




WMAP best-fit ACDM model 


1206.6 




1113.9 


Hinshaw et al. cut sky 


213.4 


2.86% 


1038.7 


WILC5YR 


242.7 


3.79% 


1053.2 


HILC5YR 


306.2 


5.75% 


1104.8 



In Table IIII1 the quadrupole and octupole powers of 
HILC5YR are shown with those of the WMAP 5 year 
best-fit ACDM model and WILC5YR. 



The angular separation between the preferred axis 
of the quadrupole anisotropy and that of the oc- 
tupole anisotropy is 2° for WILC5YR, while 12°. 1 for 
HILC5YR. The corresponding probabilities of getting 
such an alignment is 0.058% for WILC5YR, while 2.17% 
for HILC5YR. The anisotropy of HILC5YR on low mul- 
tipoles (2 < I < 5) are shown in Fig. [HI [El [16] and [17] 
with those of WILC5YR and difference maps. 




FIG. 14: Quadrupole Anisotropy [fiK]: HILC5YR (top), 
WILC5YR (middle), WILC5YR - HILC5YR (bottom) 

The linear weights of HILC5YR are shown in Fig. 
[T8l In Fig. [191 we show the variance of the lin- 
ear weights of HILC5YR, which is computed by W\ = 

(2^ + i)" 1 E„ l K m l 2 - 



VIII. COMPUTATIONAL ISSUE 

The computation of linear weights by our harmonic 
variance minimization can be split as follows: 

1) Computing Ji(h,mi,h, m 3 ) by Eq. [TB] 

2) Computing ay m , lm by Eq. M 

3) Solving the system of linear equations given by Eq. 

M 

Let's assume that there are /„ band maps of high SNR 
up to some multipole L and the cutoff multipole for linear 
weights is set to I. With the recurrence relation we 
are at present able to compute Wigner 3j symbols in Eq. 
[TBI up to high multipolcs (~ 700) fast enough. Therefore, 
step 1) put relatively little computational load. 
Step 2) requires 0(J\f 2 M) and step 3), which involves 
N by Af matrix inversion, requires 0(Af 3 ), where Af = 
(I + l) 2 /„ and M = (L + l) 2 . I is chosen to be much 
smaller than L to keep matrices in Eq. IA5I numerically 



FIG. 15: Octupole Anisotropy [/iK]: HILC5YR (top), 
WILC5YR (middle), WILC5YR - HILC5YR (bottom) 




FIG. 16: Anisotropy of I = 4 [/xK]: HILC5YR (top), 
WILC5YR (middle), WILC5YR - HILC5YR (bottom) 



non-singular (e.g. In the application to WMAP data, 
we have set I — 7 and L — 300, with consideration of 
numerical singularity and the WMAP band map's SNR.) 
Since Af = (l + l) 2 f n is much smaller than M = (i + 1) 2 , 
the total computing time is 0(Af 2 M). 

As described previously, we have found l cu toS phe- 
nomenologically by increasing it until numerical singu- 



FIG. 17: Anisotropy of I = 5: HILC5YR (top), WILC5YR 
(middle), WILC5YR - HILC5YR (bottom) 

larity in Eq. IA5I emerges. While ideally Cutoff can be 
as high as L, the optimal value of Z C utoff for the WMAP 
data seems to be much smaller than L. The numeri- 
cal singularity for (2 cu toff < L) may be attributed to 
large bandwidth and relatively small separation of the 
WMAP frequency channels, because in such configura- 
tions the frequency spectrum of foregrounds may not be 
numerically distinct enough over the channels. Through 
simple extrapolation by comparing the beamwidth and 
SNR of the Planck surveyor with those of the WMAP, 
we may make rough estimate that HILC method with 
'cutoff > 100 may be numerically stable for the Planck 
temperature data. 



IX. CONCLUSION 

In spite of the warning from the WMAP team against 
serious use of the ILC map, the ILC map has been widely 
used especially for low multipole anisotropy study, since 
the template-subtraction maps are not suitable for a 
whole sky map due to heavy foreground contamination 
within Kp2 cut. We have summarized the advantages 
and disadvantages of ILC variants and the template- 
fitting method in appendix [Cl 

For these reason, we have pursued the ILC method 
and extended it by improving the causes related to the 
complication in noise properties (e.g. smoothing on the 
disjoint regions, Monte-Carlo 'bias' correction). Through 
our effort toward improved ILC implementation, we have 
developed a harmonic variance minimization method, 
which is derived by converting a general form of pixel- 
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FIG. 18: The HILC5YR linear weights for the K, Ka, Q, V, 
and W band map (from top to bottom) 

domain approach into spherical harmonic space. In our 
approach, spatial variability of linear weights is incor- 
porated in a self-contained manner and linear weights 
are continuous over whole sky. Thanks to full corre- 
spondence to a general pixel-domain method, physical 
interpretation on our method is quite straightforward. 
In variance minimization, there exists a cross term be- 
tween residual foreground and CMB, which makes the 
linear combination of minimum variance differ from that 
of minimum foreground. We have developed an iterative 
method, where perturbative correction is made for the 
cross term. 

The simulations showed that our method yields reli- 
able and stable reconstruction of the CMB also on low- 
est multipoles. By applying it to the WMAP data, 
we have obtained a CMB map, whose power spectra 
makes better agreement with the WMAP best-fit A CDM 
model than the WILC5YR. The CMB map and lin- 
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FIG. 19: Variance of HILC5YR linear weights 

ear weights, which we have obtained, are available from 
http : //www.nbi . dk/~jkim/hilc. 

Capacity of our method to clean foreground, whose 
spectral indice is a rapidly varying function of position, 
is limited by the SNR of data on high multipoles, since 
increasing the cutoff multipole of linear weights requires 
more data to avoid numerical singularity in Eq. IA5I 
When the data of good SNR on high multipoles are avail- 
able, we will be able to make better reduction of fore- 
grounds by setting the cutoff multipole to higher value. 
This is similar to saying that we are able to assume finer 
disjoint regions in pixel-domain approach, when data of 
better pixel resolution are available. 

The foreground reduction method by template fitting 
is unable to take full advantage of the high resolution 
low noise Planck data, since the available templates do 
not have such high resolution and low noise over a whole 
sky. Unlike the template method, the effectiveness of 
our method scales with frequency channel numbers and 
angular resolution of the observation data. Hence our 
method is more suitable for the Planck surveyor, which 
has nine frequency channels with good SNR and angular 
resolution. 

Our method, which is presented for the application to 
temperature data, may be easily extended to the polar- 
ization data by making the following replacement in Eq. 

us 

(h h h\ (h h h \ 

^000,! ^ ^ ±2 =p2 J ' 

We believe a blind approach like our method is desirable 
for the analysis of polarization data from the upcoming 
Planck satellite, since the availability of polarized fore- 
ground templates is quite limited. For these reasons, the 
method presented in this paper will stand out among sev- 
eral foreground reduction methods, when low-noise high 
resolution data of CMB temperature and polarization are 
available from the Planck surveyor (l7j . 
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APPENDIX A: COMPUTING LINEAR WEIGHTS 
VIA MATRIX OPERATIONS 



APPENDIX B: THE EFFECT OF THE CROSS 
TERM IN VARIANCE MINIMIZATION 

To better understand the effect of the cross term, we 
investigate a simple case of two frequency observation 
channels and foregrounds of uniform frequency spectra. 
Since the analysis in Sec. [V] and I VII are carried out on 
the multipolc of high SNR, we neglect instrument noise. 
We assume linear weights to be constant, because the 
foreground spectra are assumed to be spatially uniform. 
Keeping the CMB signal unchanged, we assign constant 
linear weights w and (1 — w) for the frequency channel 1 
and 2 respectively. Since spherical harmonic coefficients 
of the linear combination map is 



We present the solutions of Eq. ITTT [T2l and [T4l through 
matrix operations. Consider n linear weights for maps of 
k frequency channels. In matrix notation, the constraints 
given by Eq. [TT1 and [L2l are as follows: 



n • w = e. 



(Al) 



ah 



wa\t + {l-w)aZ+a<i 



fg2 



_cmb 



the variance of the linear combination map is as follows: 



II is a r x n matrix, given by 



n, 



1 : k(i - 1) + 1 < j < ki 
: otherwise 



and e is a column vector of length n/k, given by 

V4n : j = 1 
: j>l 

In matrix notation, Eq. [T4"lis as follows: 



A w = II L, 



(A2) 



where 



l l'm'lmi 



W, = W 



Lj» = Xl'm'i 

for j' = k(l' 2 + I' + m') + i', j = k(l 2 + l + m)+i, and 
j" = I' 2 + I' + m! . A is a n x n matrix, and w and L are 
column vectors of length n, and length n/k respectively. 
With Eq. IA2i w is solved in terms of n/k undetermined 
Langrange multipliers, provided that A is invertible: 



A _1 n T L 



(A3) 



With Eq. lAll and |A31 the undetermined n/k Langrange 
multipliers are given by: 



l = (nA x n 



1 ttTn-1 



e, 



(A4) 



Therefore, linear weights in spherical harmonic space is 
given as follows: 



w = A- 1 n T (nA- i n T )- i e. 



-1 ttT\-1. 



(A5) 



Eq. IA5I is not reduced to w = II x e, since II and n T 
are not square matrices. 



where off* / and tiiff,/ denote the spherical harmonic co- 
efficients of foregrounds at frequency channel 1 and 2 
respectively. Since the value of w, which minimizes Eq. 
iBil is 



-Eae[(aft-^)(*+* 

I'm' L 



V 



I'm' 



a v 



fe2 ,2 



the linear combination map of minimum variance has the 
following spherical harmonic coefficients: 



^cmb 



1 



fg2 



-(a 



a zf ) J2 Re 



(a 



fei 




The power spectra of our galactic foregrounds has uneven 
distribution with high concentration on lowest multipole 
and azimuthal mode (i.e. a 2 g 3> o,<f m i)- Therefore, the 
spherical harmonic coefficients of the linear combination 
map is 



air. 



nmb 



fg2 

u lm ) f K l 



"20 



fg2 
a 2Q 



(B2) 



According to Eq. [HI a lm w for (I = 2, to = 0). This 
is in agreement with the suppression on lowest multi- 
pole powers in an Internal Linear Combination (ILC) 
map, reported by [§, [lj|. Therefore, we attribute the 
existence of the cross term and highly uneven power 
spectrum of foregrounds to the suppression on low mul- 
tipole power in an Internal Linear Combination (ILC) 
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map. To test our hypothesis, we have carried out 
simulations with foregrounds of flat power spectrum, 
which are derived from the WMAP Maximum-Entropy- 
Method(MEM) foregrounds as follows: 



2 fg ' 

l lm 



la U l 



(B3) 



where a£ is spherical harmonic coefficients of the 
WMAP MEM foreground and Cf s = £ m \af m \ 2 /{2l + 1). 



The role of y E ; Cf s /('i — Zq + 1) is to match the total 
power of the flat foreground with that of the MEM fore- 
grounds. The result of simulations show that there is 
no observable suppression on low multipole powers when 
foregrounds have flat power spectrum, thereby support- 
ing our hypothesis. 

This suppression on low multipole power can be re- 
duced by excluding the multipoles of high foreground 
concentration (i.e. low multipoles) in variance minimiza- 
tion process. Consider minimizing the following variance 
where low multipoles are excluded: 



E 

l'>lo 



fgi 



^cmb 



where foreground power in the multipole range (I' < Iq) 
are much bigger than those on higher multipoles. Then, 
the linear combination map has the following spherical 
harmonic coefficients: 



a/7 



cmb 
a lm 



1 



E 

1'>Io,itl' 



fgl 

fym' 



fg2 
l Vm> 



(B4) 



fg2 



E 



fg2 



(a 



fgi 



,fg2\ 



E 

1'>Iq ,m' 



Re 



(a fgl - a fg2 1 (a ts2 

\ a l>m> a l'm>) y a l'm> 



)* 



With enough number of summation terms, the term in- 
side the big parenthesis of Eq. IB4I approaches zero be- 
cause of cancellation. If we had not excluded low multi- 
poles, the cancellation would be ineffective, due to asym- 
metric distribution of foreground power (i.e. high con- 
centration on low multipoles.). Eq. IB4I is valid for a im 
(I < lo) as well as (I > Iq). In other words, though the 
linear weights were determined through variance mini- 
mization over (Iq < I'), they are applicable to foreground 
reduction in ai m (I < lo). 

We have also investigated the non-uniform frequency 
spectra case by resorting to a numerical investigation. 
We find that there is a tendency of suppression over 
multipole range (l - l cuto s < I < l + Cutoff), where 



'cutoff is the assumed cutoff multipole of w\ m . We also 
find that the suppression is reduced by excluding the 
multipoles of high foreground concentration in variance 
minimization process, just as in the uniform frequency 
spectra case. 



APPENDIX C: COMPARISON OF THE ILC 
METHOD VARIANTS AND THE 
TEMPLATE-FITTING METHOD 

In this section, we are discussing briefly the advan- 
tages and disadvantages of ILC variants and a template 
fitting method. The template fitting method is the fore- 
ground reduction method most importantly employed by 
the WMAP team. It has advantage that it has less 
complicated noise properties and is free from the Cos- 
mic Covariance problem, while it has the disadvantages 
that it relies on foreground templates of external sources, 
hence requiring extrapolation to observation frequencies 
and currently unable to provide a whole sky map, due to 
heavy foreground contamination within Kp2 cut. WILC 
is the ILC implementation by the WMAP team. It has 
the advantages that it utilizes the boundary shape infor- 
mation of galactic foregrounds and scales with the num- 
ber of observation frequency channels, while its disadvan- 
tages are sharp boundaries of disjoint regions (Smooth- 
ing boundaries in the final map making does not solve 
discontinuity problem completely, since region definition 
with sharp boundary are used in variance minimization.), 
dubious bias correction of the 'Cosmic Covariance' by 
Monte-Carlo CMB, its reliance on the pre-defined dis- 
joint regions (not being a completely blind approach). 
TCM3YR Q is the ILC variant, where the variance of 
each multipole is minimized separately. It has the ad- 
vantages that the dependency of foreground power on 
angular scales is reflected and scales with the number 
of frequency channels, while it has the disadvantages of 
sharp boundaries, 'Cosmic Covariance' problem, need for 
the pre-defined disjoint regions. SILC3YR [18] is another 
ILC variant. Instead of using disjoint region definition by 
the WMAP team, disjoint regions were derived from the 
MEM reconstructed foregrounds. It has the advantages 
that the definition of the disjoint regions may be optimal 
than those of WILC, and it scales with the number of fre- 
quency channels, while it has the disadvantages of sharp 
boundaries, 'Cosmic Covariance' problem, need for the 
MEM foreground (traces MEM foregrounds and hence 
WILC3YR). HILC, which is our ILC implementation, has 
the advantages that it does not rely on the definition of 
disjoint regions (hence no sharp boundaries), scalability 
with the number of observation frequency channels and 
angular resolution of data, while it has the disadvantages 
that it is computationally intensive, and has the 'Cosmic 
Covariance' problem, though reduced. 
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